###################################################################
### Political Institutions, Energy Transitions, and Air Quality ### 
### Evidence from Global Urban Areas                            ###
### Main empirical findings / difference in differences         ###                  
### This version: November 16, 2024                             ###                   
###################################################################

#The code below reproduces the descriptive figures in the main text and SI

#---------------------------------------------------#
# Preliminaries                                     #
#---------------------------------------------------#

rm(list = ls())

#List of packages for session
.packages = c("tidyverse", "did", "texreg", "countrycode",
              "fixest", "WDI", "gridExtra", "sf", "nngeo")

#Install CRAN packages (if not already installed)
.inst = .packages %in% installed.packages()
if(length(.packages[!.inst]) > 0) install.packages(.packages[!.inst])

#Loading packages into session 
lapply(.packages, require, character.only=TRUE)

# Set Working Directory to wherever files are downloaded
#setwd("")

#---------------------------------------------------#
# Part 1: Loading the data                          #
#---------------------------------------------------#

master_cities_panel <- read.csv("master_cities_panel.csv")
master_pollution_dem <- read.csv("master_pollution_dem.csv")
master_cities_panel_robustness <- read.csv("master_cities_panel_robustness.csv")

#Coal power plants and air pollution original raw data
coal_plants <- read.csv("coal_power_plants.csv")
air_pollution_global2 <- read.csv("air_pollution_global2.csv")

##############################
###### Figures Main Text #####
##############################

#---------------------------------------------------#
# Figure 1: Identification strategy                 #
#---------------------------------------------------#

##1.2 Checking duplicates (the original dataset has some duplicates, i.e. Rochester NY)

dup <- air_pollution_global2 %>%
  group_by(ID, City, Country, Year) %>%
  summarise(tot = n()) %>%
  filter(tot > 1) #several cities in the United States
#Salem, Richmond, Lakewood repeated 1; Rochester (2) and Bloomington (2) repeated 2
# this is equal to 3*20 + 4*40 = 220 observations repeated for a total of 260020-259800

##1.3. Eliminating duplicates

dta <- air_pollution_global2 %>%
  distinct(ID, City, Country, Year, .keep_all = T) %>% #259,800 unique observations
  mutate(country2 = countrycode(Country, 'country.name', 'iso3c'),
         latitude = Latitude,
         longitude = Longitude)

##1.4 Transform the coal power plant data

coal_plants <- st_as_sf(coal_plants, coords = c("Longitude", "Latitude"), crs = 4326)

##1.5 Geo-referencing the city data using the coal power plants coordinates

dta_geo <- dta %>% 
  mutate(city_uid = paste(City, country2, sep = "-")) %>%
  dplyr::select(ID, city_uid, country2, Latitude, Longitude) %>% distinct(ID, city_uid, country2, Latitude, Longitude, .keep_all = F)

dta_geo2 <- master_pollution_dem %>%
  distinct(ID, border_city_t1, border_city_t2, .keep_all = F) %>%
  left_join(dta_geo, by = "ID") %>%
  mutate(border_city_t3 = ifelse(border_city_t1 == 1 & border_city_t2 == 1, 1, 0)) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = st_crs(coal_plants)) 

#1.6 Create the plot

world_map <- map_data("world")

pdf(file = "figures/cities_dem.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 6) # The height of the plot in inches
ggplot() +
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = subset(dta_geo2, border_city_t1 == 1),
    colour = "#CD27D8",
    shape = 17,
    size = 0.65,
    alpha = 0.5
  ) +
  geom_sf(
    data = subset(dta_geo2, border_city_t2 == 1),
    colour = "darkgreen",
    shape = 17,
    size = 0.65,
    alpha = 0.5,
  ) +
  geom_sf(
    data =subset(dta_geo2, border_city_t3 == 1),
    colour = "#278AD8",
    shape = 17,
    size = 0.65,
    alpha = 0.5,
  ) +
  theme(legend.position="none")
dev.off()

#---------------------------------------------------#
# Figure 2: Exposure to coal-power infrastructure.  #
#---------------------------------------------------#

##2.1: Separating coal power plants by their status 

coal_plants_operating <- coal_plants %>%
  filter(Status %in%c("operating","construction")) %>%
  mutate(unit = row_number()) 

coal_plants_retired <- coal_plants %>%
  filter(Status %in% c("retired","mothballed","shelved")) %>%
  mutate(unit = row_number()) 

#Here I am separating coal power plants by overall capacity

coal_plants_retired_high <- coal_plants %>%
  filter(Status %in% c("retired","mothballed","shelved")) %>%
  filter(Capacity..MW. > 300.0) %>%
  mutate(unit = row_number()) 

coal_plants_retired_low <- coal_plants %>%
  filter(Status %in% c("retired","mothballed","shelved")) %>%
  filter(Capacity..MW. < 300.0) %>%
  mutate(unit = row_number()) 

coal_plants_counterfactual <- coal_plants %>% #It refers to plants that were never constructed
  filter(Status %in% c("permitted","announced","cancelled")) %>%
  mutate(unit = row_number()) 

##2.2 Cross-section of cities

cross_section_cities <- dta_geo2 %>%
  distinct(ID, city_uid, country2, geometry, .keep_all = FALSE)

##2.3 Calculating distance to power plants

###Counterfactual

cities_coal_plants_counterfactual <- st_nn(cross_section_cities, coal_plants_counterfactual, k = 1, returnDist = T)  %>% 
  set_names("unit", "distance_km_count") %>% 
  map_df(., unlist)

cities_coal_plants_counterfactual2 <- cbind(cross_section_cities, cities_coal_plants_counterfactual) %>%
  as.data.frame() %>%
  left_join(coal_plants_counterfactual, by = "unit") 

###Operating

cities_coal_plants_operating <- st_nn(cross_section_cities, coal_plants_operating, k = 1, returnDist = T)  %>% 
  set_names("unit", "distance_km_oper") %>% 
  map_df(., unlist)

cities_coal_plants_operating2 <- cbind(cross_section_cities, cities_coal_plants_operating) %>%
  as.data.frame() %>%
  left_join(coal_plants_counterfactual, by = "unit") 

###Retired

cities_coal_plants_retired <- st_nn(cross_section_cities, coal_plants_retired, k = 1, returnDist = T)  %>% 
  set_names("unit", "distance_km_ret") %>% 
  map_df(., unlist)

cities_coal_plants_retired2 <- cbind(cross_section_cities, cities_coal_plants_retired) %>%
  as.data.frame() %>%
  left_join(coal_plants_retired, by = "unit") 

###Retired high

cities_coal_plants_retired_high <- st_nn(cross_section_cities, coal_plants_retired_high, k = 1, returnDist = T)  %>% 
  set_names("unit", "distance_km_ret") %>% 
  map_df(., unlist)

cities_coal_plants_retired_high2 <- cbind(cross_section_cities, cities_coal_plants_retired_high) %>%
  as.data.frame() %>%
  left_join(coal_plants_retired_high, by = "unit")

###Retired low

cities_coal_plants_retired_low <- st_nn(cross_section_cities, coal_plants_retired_low, k = 1, returnDist = T)  %>% 
  set_names("unit", "distance_km_ret") %>% 
  map_df(., unlist)

cities_coal_plants_retired_low2 <- cbind(cross_section_cities, cities_coal_plants_retired_low) %>%
  as.data.frame() %>%
  left_join(coal_plants_retired_low, by = "unit") 

##2.4: Merging the distance calculations

###50k

cities_close_retired50k <- cities_coal_plants_retired2 %>%
  filter(distance_km_ret < 50000) %>% 
  mutate(close_retired50k = 1) %>% 
  rename(distance_km_ret50k = distance_km_ret,
         Retired.year50k = Retired.year) %>%
  dplyr::select(ID, close_retired50k, distance_km_ret50k, Retired.year50k)

cities_close_retired50k_high <- cities_coal_plants_retired_high2 %>%
  filter(distance_km_ret < 50000) %>% 
  mutate(close_retired50k_high = 1) %>% 
  rename(distance_km_ret50k_high = distance_km_ret,
         Retired.year50k_high = Retired.year) %>%
  dplyr::select(ID, close_retired50k_high, distance_km_ret50k_high, Retired.year50k_high)

cities_close_retired50k_low <- cities_coal_plants_retired_low2 %>%
  filter(distance_km_ret < 50000) %>% 
  mutate(close_retired50k_low = 1) %>% 
  rename(distance_km_ret50k_low = distance_km_ret,
         Retired.year50k_low = Retired.year) %>%
  dplyr::select(ID, close_retired50k_low, distance_km_ret50k_low, Retired.year50k_low)

cities_close_operating50k <- cities_coal_plants_operating2 %>%
  filter(distance_km_oper < 50000) %>% 
  mutate(close_operating50k = 1) %>% 
  rename(distance_km_oper50k = distance_km_oper) %>%
  dplyr::select(ID, close_operating50k, distance_km_oper50k)

cities_close_counterfactual50k <- cities_coal_plants_counterfactual2 %>%
  filter(distance_km_count < 50000) %>% 
  mutate(close_counterfactual50k = 1) %>% 
  rename(distance_km_count50k = distance_km_count) %>%
  dplyr::select(ID, close_counterfactual50k, distance_km_count50k)

###100k

cities_close_retired100k <- cities_coal_plants_retired2 %>%
  filter(distance_km_ret < 100000) %>% 
  mutate(close_retired100k = 1) %>% 
  rename(distance_km_ret100k = distance_km_ret,
         Retired.year100k = Retired.year) %>%
  dplyr::select(ID, close_retired100k, distance_km_ret100k, Retired.year100k)

cities_close_operating100k <- cities_coal_plants_operating2 %>%
  filter(distance_km_oper < 100000) %>% 
  mutate(close_operating100k = 1) %>% 
  rename(distance_km_oper100k = distance_km_oper) %>%
  dplyr::select(ID, close_operating100k, distance_km_oper100k)

cities_close_counterfactual100k <- cities_coal_plants_counterfactual2 %>%
  filter(distance_km_count < 100000) %>% 
  mutate(close_counterfactual100k = 1) %>% 
  rename(distance_km_count100k = distance_km_count) %>%
  dplyr::select(ID, close_counterfactual100k, distance_km_count100k)

###150k

cities_close_retired150k <- cities_coal_plants_retired2 %>%
  filter(distance_km_ret < 150000) %>% 
  mutate(close_retired150k = 1) %>% 
  rename(distance_km_ret150k = distance_km_ret,
         Retired.year150k = Retired.year) %>%
  dplyr::select(ID, close_retired150k, distance_km_ret150k, Retired.year150k)

cities_close_operating150k <- cities_coal_plants_operating2 %>%
  filter(distance_km_oper < 150000) %>% 
  mutate(close_operating150k = 1) %>% 
  rename(distance_km_oper150k = distance_km_oper) %>%
  dplyr::select(ID, close_operating150k, distance_km_oper150k)

cities_close_counterfactual150k <- cities_coal_plants_counterfactual2 %>%
  filter(distance_km_count < 150000) %>% 
  mutate(close_counterfactual150k = 1) %>% 
  rename(distance_km_count150k = distance_km_count) %>%
  dplyr::select(ID, close_counterfactual150k, distance_km_count150k)

###450k

cities_close_retired450k <- cities_coal_plants_retired2 %>%
  filter(distance_km_ret < 450000) %>% 
  mutate(close_retired450k = 1) %>% 
  rename(distance_km_ret450k = distance_km_ret,
         Retired.year450k = Retired.year) %>%
  dplyr::select(ID, close_retired450k, distance_km_ret450k, Retired.year450k)

cities_close_operating450k <- cities_coal_plants_operating2 %>%
  filter(distance_km_oper < 450000) %>% 
  mutate(close_operating450k = 1) %>% 
  rename(distance_km_oper450k = distance_km_oper) %>%
  dplyr::select(ID, close_operating450k, distance_km_oper450k)

cities_close_counterfactual450k <- cities_coal_plants_counterfactual2 %>%
  filter(distance_km_count < 450000) %>% 
  mutate(close_counterfactual450k = 1) %>% 
  rename(distance_km_count450k = distance_km_count) %>%
  dplyr::select(ID, close_counterfactual450k, distance_km_count450k)

#Double-check distance worked

cities_geo_distance <- cross_section_cities %>%
  left_join(cities_close_operating50k, by = "ID") %>%
  left_join(cities_close_retired50k, by = "ID") %>%
  left_join(cities_close_retired50k_high, by = "ID") %>%
  left_join(cities_close_retired50k_low, by = "ID") %>%
  left_join(cities_close_operating100k, by = "ID") %>%
  left_join(cities_close_retired100k, by = "ID") %>%
  left_join(cities_close_operating150k, by = "ID") %>%
  left_join(cities_close_retired150k, by = "ID") %>%
  left_join(cities_close_operating450k, by = "ID") %>%
  left_join(cities_close_retired450k, by = "ID") %>%
  filter(!ID %in% c("3168")) #There is an error with the coordinates of this city in the original data

#Left-Panel

# Define the latitude and longitude extent
xmin <- -25
xmax <- 40
ymin <- 30
ymax <- 75

pdf(file = "figures/europe.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 12) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("France","Belgium","Germany","Switzerland",
                                                      "Netherlands","United Kingdom", "Italy", "Spain",
                                                      "Denmark","Sweden","Norway","Austria","Poland",
                                                      "Ukraine","Latvia","Estonia","Lithuania",
                                                      "Hungary","Croatia","Serbia","Greece","Albania",
                                                      "Romania","Czech Republic","Slovakia",
                                                      "Bosnia and Herzegovina","Slovenia", "Portugal",
                                                      "Bulgaria","Belarus","Kosovo","North Macedonia",
                                                      "Moldova","Finland","Turkey","UK","Ireland",
                                                      "Morocco","Russia","Israel","Syria",
                                                      "Libya","Algeria","Tunisia","Egypt")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#E11ECB",
    size = 1.2,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, is.na(close_operating50k)),
    colour = "#CBE11E",
    size = 0.75,
  ) +
  geom_sf(
    data = subset(coal_plants, Status == "operating"),
    colour = "#1ECBE1",
    shape = 17,
    size = 1.2,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) + 
  theme(legend.position="none")
dev.off()

#Right-Panel

# Define the latitude and longitude extent
xmin <- 75
xmax <- 130
ymin <- 15
ymax <- 55

pdf(file = "figures/asia.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 12) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("China","India",
                                                      "Nepal","Myanmar",
                                                      "Laos","Thailand",
                                                      "Vietnam","Japan",
                                                      "Mongolia","Russia",
                                                      "Kazakhstan","South Korea",
                                                      "North Korea")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#E11ECB",
    size = 1.2,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, is.na(close_operating50k)),
    colour = "#CBE11E",
    size = 0.75,
  ) +
  geom_sf(
    data = subset(coal_plants, Status == "operating"),
    colour = "#1ECBE1",
    shape = 17,
    size = 1.2,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) + 
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure 3: Temporal trends in mortality rates across countries by income category.  #
#------------------------------------------------------------------------------------#

gdppc_atlas <- WDI(country = "all", indicator = "NY.GNP.PCAP.CD", start = 2000,
                   end = 2019, extra = FALSE, cache = NULL)

gdppc_atlas <- gdppc_atlas %>%
  filter(!country %in% c("Arab World","Early-demographic dividend","East Asia & Pacific",
                         "East Asia & Pacific (excluding high income)",
                         "East Asia & Pacific (IDA & IBRD countries)",
                         "Europe & Central Asia (excluding high income)",
                         "Europe & Central Asia (IDA & IBRD countries)",
                         "Fragile and conflict affected situations",
                         "Heavily indebted poor countries (HIPC)", "High income",
                         "IBRD only","IDA & IBRD total","IDA blend","IDA only",
                         "IDA total","Late-demographic dividend",
                         "Latin America & Caribbean",
                         "Latin America & Caribbean (excluding high income)",
                         "Latin America & the Caribbean (IDA & IBRD countries)",
                         "Least developed countries: UN classification",
                         "Low & middle income","Low income","Lower middle income",
                         "Middle East & North Africa",
                         "Middle East & North Africa (excluding high income)",
                         "Middle East & North Africa (IDA & IBRD countries)",
                         "Middle income","North America","Not classified",
                         "OECD members","Other small states","Pacific island small states",
                         "Post-demographic dividend","Pre-demographic dividend",
                         "Small states","South Asia","South Asia (IDA & IBRD)",
                         "Sub-Saharan Africa","Sub-Saharan Africa (excluding high income)",
                         "Sub-Saharan Africa (IDA & IBRD countries)", "Upper middle income",
                         "World", "Caribbean small states", "European Union",
                         "Europe & Central Asia", "Euro area",
                         "Central Europe and the Baltics",
                         "West Bank and Gaza","Virgin Islands (U.S.)",
                         "Vanuatu","Tuvalu","Turks and Caicos Islands","Tonga",
                         "St. Vincent and the Grenadines","St. Martin (French part)",
                         "St. Lucia","St. Kitts and Nevis","Solomon Islands",
                         "Sint Maarten (Dutch part)","Palau","Northern Mariana Islands",
                         "New Caledonia","Nauru","Micronesia, Fed. Sts.",
                         "Marshall Islands","Macao SAR, China","Kiribati",
                         "Isle of Man","Hong Kong SAR, China","Greenland","Gibraltar",
                         "French Polynesia","Fiji","Faroe Islands","Dominica",
                         "Channel Islands","Cayman Islands","British Virgin Islands",
                         "Bermuda","Antigua and Barbuda","American Samoa",
                         "Grenada","Aruba","Curacao")) %>%
  mutate(country2 = countrycode(iso2c, 'iso2c', 'iso3c')) %>%
  group_by(country2) %>%
  summarise(gdppc_atlas = mean(NY.GNP.PCAP.CD, na.rm = T)) %>%
  filter(!is.na(country2) & !is.nan(gdppc_atlas))

income_cat <- gdppc_atlas %>%
  group_by(country2) %>%
  summarise(gdppc_atlas = mean(gdppc_atlas, na.rm = T)) %>%
  mutate(income_cat = case_when(gdppc_atlas <= 1145 ~ "Low Income",
                                gdppc_atlas > 1145 & gdppc_atlas <= 4515 ~ "Low-Middle Income",
                                gdppc_atlas > 4515 & gdppc_atlas <= 14005 ~ "Upper-Middle Income",
                                gdppc_atlas > 14005 ~ "High Income",
                                TRUE ~ NA_character_))

air_pollution_global_final <- dta %>%
  left_join(income_cat, by = "country2") #Missing countries for income, DPRK, Curacao, Fiji, French P, New Caledonia, Palesine, Jersey, French Guiana, Tauwan, Guadalpue, Martinique, Salomon Islands

#Summarizing the data

plot1_dat <- air_pollution_global_final %>%
  filter(!is.na(income_cat)) %>%
  group_by(income_cat, Year) %>%
  summarise(Rates_PM = mean(Rates_PM, na.rm = T),
            Rates_NO2 = mean(Rates_NO2, na.rm = T),
            Rates_O3 = mean(Rates_O3, na.rm = T)) %>%
  ungroup() 

plot1a <- ggplot(plot1_dat, aes(x = Year, y = Rates_PM, group = income_cat, color = income_cat)) +
  geom_line() +
  labs(title = "Average Mortality Associated with PM2.5",
       x = "Year",
       y = "Cases Attributable per 100,000",
       color = "Income Group") +
  scale_color_brewer(palette = "Set2") +
  theme_classic()

plot1b <- ggplot(plot1_dat, aes(x = Year, y = Rates_NO2, group = income_cat, color = income_cat)) +
  geom_line() +
  labs(title = "Average Mortality Associated with Nitrogen Dioxide",
       x = "Year",
       y = "Cases Attributable per 100,000",
       color = "Income Group") +
  scale_color_brewer(palette = "Set2") +
  theme_classic()

pdf(file = "figures/fig1_incomegroups.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot1a, plot1b, nrow = 2)  
dev.off()

#-----------------------------------------------------------------------------------------------------------------------------#
# Figure 4: Countries with the largest absolute increases and decreases in cases attributable to PM2.5 and nitrogen dioxide.  #
#-----------------------------------------------------------------------------------------------------------------------------#

dta_change <- air_pollution_global_final %>%
  group_by(City, ID, country2, Latitude, Longitude) %>%
  filter(Year %in% c(2000, 2019)) %>%
  mutate(change_pm_rate = Rates_PM-lag(Rates_PM),
         change_no2_rate = Rates_NO2-lag(Rates_NO2),
         change_o3_rate = Rates_O3-lag(Rates_O3),
         change_co2 = CO2-lag(CO2)) %>%
  filter(!is.na(change_pm_rate)) %>%
  dplyr::select(City, ID, country2, Rates_PM, Rates_O3, Rates_NO2, CO2, 
                change_pm_rate, change_no2_rate, change_o3_rate, change_co2)

plot2_dat <- dta_change %>%
  group_by(country2) %>%
  summarise(change_pm_rate = mean(change_pm_rate, na.rm = T),
            change_no2_rate = mean(change_no2_rate, na.rm = T),
            change_o3_rate = mean(change_o3_rate, na.rm = T),
            change_co2 = mean(change_co2, na.rm = T)) %>%
  filter(country2 != "PRI")

top_10_pm25 <- plot2_dat %>%
  arrange(desc(change_pm_rate)) %>%          
  slice_head(n = 10) 

top_10_no2 <- plot2_dat %>%
  arrange(desc(change_no2_rate)) %>%          
  slice_head(n = 10) 

bottom_10_pm25 <- plot2_dat %>%
  arrange(change_pm_rate) %>%          
  slice_head(n = 10) 

bottom_10_no2 <- plot2_dat %>%
  arrange(change_no2_rate) %>%          
  slice_head(n = 10) 

plot2a <- ggplot(top_10_pm25, aes(x = reorder(country2, -change_pm_rate), y = change_pm_rate)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Countries with largest increase in PM2.5",
       x = "Country", y = "Avg Pct Ch 00-19") +
  theme_classic()

plot2b <- ggplot(top_10_no2, aes(x = reorder(country2, -change_no2_rate), y = change_no2_rate)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Countries with largest increase in NO2",
       x = "Country", y = "Avg Pct Ch 00-19") +
  theme_classic()

plot2c <- ggplot(bottom_10_pm25, aes(x = reorder(country2, -change_pm_rate), y = change_pm_rate)) +
  geom_bar(stat = "identity", fill = "salmon") +
  labs(title = "Countries with largest decrease in PM2.5",
       x = "Country", y = "Avg Pct Ch 00-19") +
  theme_classic()

plot2d <- ggplot(bottom_10_no2, aes(x = reorder(country2, -change_no2_rate), y = change_no2_rate)) +
  geom_bar(stat = "identity", fill = "salmon") +
  labs(title = "Countries with largest decrease in NO2",
       x = "Country", y = "Avg Pct Ch 00-19") +
  theme_classic()

pdf(file = "figures/fig2_top_changes.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot2a, plot2b, plot2c, plot2d, nrow = 2)  
dev.off()

#------------------------------------------------------------------------------------------------------#
# Figure 5: Temporal trends in mortality rates from PM2.5 and nitrogen dioxide in selected countries.  #
#------------------------------------------------------------------------------------------------------#

timeseries_dat <- air_pollution_global_final %>%
  group_by(Country, Year) %>%
  summarise(avgpm = mean(Rates_PM, na.rm = T),
            avgno2 = mean(Rates_NO2, na.rm = T),
            avgo3 = mean(Rates_O3, na.rm = T)) %>%
  filter(Country %in% c("United States","Mexico",
                        "Russia","India","China", "Brazil",
                        "Nigeria","France","Japan"))

plot4a <- ggplot(timeseries_dat, aes(x = Year, y = avgpm, group = Country, color = Country)) +
  geom_line() +
  labs(title = "Average Mortality Associated with PM2.5",
       x = "Year",
       y = "Average Cases Attributable per 100,000",
       color = "Country") +
  scale_color_brewer(palette = "Set2") +
  theme_classic()

plot4b <- ggplot(timeseries_dat, aes(x = Year, y = avgno2, group = Country, color = Country)) +
  geom_line() +
  labs(title = "Average Mortality Associated with NO2",
       x = "Year",
       y = "Average Cases Attributable per 100,000",
       color = "Country") +
  scale_color_brewer(palette = "Set2") +
  theme_classic()

pdf(file = "figures/fig4_selectedcountries.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot4a, plot4b, nrow = 2)  
dev.off()

#------------------------------------------------------------------------------------#
# Figure 6: Cross-city variation in mortality rates from air pollution.              #
#------------------------------------------------------------------------------------#

#Show variation by country USA

usa_cities <- air_pollution_global_final %>%
  filter(country2 == "USA") %>%
  group_by(City, ID, country2, Latitude, Longitude) %>%
  filter(Year %in% c(2000, 2019)) %>%
  mutate(change_pm_rate = Rates_PM-lag(Rates_PM),
         change_no2_rate = Rates_NO2-lag(Rates_NO2),
         change_o3_rate = Rates_O3-lag(Rates_O3),
         change_co2 = CO2-lag(CO2)) %>%
  filter(!is.na(change_pm_rate)) %>%
  dplyr::select(City, ID, country2, Rates_PM, Rates_O3, Rates_NO2, CO2, 
                change_pm_rate, change_no2_rate, change_o3_rate, change_co2) %>% ungroup()

top_10usa_cities <- usa_cities %>%
  arrange(desc(change_pm_rate)) %>%  # Arrange in descending order
  slice(1:10)                       # Select top 10 rows

# Select bottom 10 rows with the lowest values
bottom_10usa_cities <- usa_cities %>%
  arrange(change_pm_rate) %>%        # Arrange in ascending order
  slice(1:10)                       # Select bottom 10 rows

# Combine top and bottom 10 rows
combined_usa <- bind_rows(top_10usa_cities, bottom_10usa_cities)

#Show variation by country China

china_cities <- air_pollution_global_final %>%
  filter(country2 == "CHN") %>%
  group_by(City, ID, country2, Latitude, Longitude) %>%
  filter(Year %in% c(2000, 2019)) %>%
  mutate(change_pm_rate = Rates_PM-lag(Rates_PM),
         change_no2_rate = Rates_NO2-lag(Rates_NO2),
         change_o3_rate = Rates_O3-lag(Rates_O3),
         change_co2 = CO2-lag(CO2)) %>%
  filter(!is.na(change_pm_rate)) %>%
  dplyr::select(City, ID, country2, Rates_PM, Rates_O3, Rates_NO2, CO2, 
                change_pm_rate, change_no2_rate, change_o3_rate, change_co2) %>% ungroup()

top_10china_cities <- china_cities %>%
  arrange(desc(change_pm_rate)) %>%  # Arrange in descending order
  slice(1:10)                       # Select top 10 rows

# Select bottom 10 rows with the lowest values
bottom_10china_cities <- china_cities %>%
  arrange(change_pm_rate) %>%        # Arrange in ascending order
  slice(1:10)                       # Select bottom 10 rows

# Combine top and bottom 10 rows
combined_china <- bind_rows(top_10china_cities, bottom_10china_cities)

plot9a <- ggplot(combined_usa, aes(x = reorder(City, -change_pm_rate), y = change_pm_rate, fill = change_pm_rate > 0)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("steelblue", "salmon"), 
                    labels = c(paste("Decrease"), paste("Increase"))) +
  labs(title = "Top and Bottom 10 US Cities by Pct Change in PM2.5",
       x = "City", y = "Change", fill = "Change") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
  

plot9b <- ggplot(combined_china, aes(x = reorder(City, -change_pm_rate), y = change_pm_rate, fill = change_pm_rate > 0)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("steelblue", "salmon"), 
                    labels = c(paste("Decrease"), paste("Increase"))) +
  labs(title = "Top and Bottom 10 Chinese Cities by Pct Change in PM2.5",
       x = "City", y = "Change", fill = "Change") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

pdf(file = "figures/fig9_cityvariation.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot9a, plot9b, nrow = 2)  
dev.off()

#------------------------------------------------------------------------------------#
# Figure 7: Air pollution mortality in the 20 most populated cities in the world.    #
#------------------------------------------------------------------------------------#

global_cities <- air_pollution_global_final %>%
  filter(Year == 2019) %>%
  filter(Population > 14497914) 

global_cities2 <- air_pollution_global_final %>% #filter by city ID
  filter(ID %in% c(10745, 12141, 12456, 3915, 11920, 12911, 13017, 155, 12882,
                   12445, 3686, 10774, 3574, 14, 950, 9933, 1308, 6887, 6998, 9752)) %>%
  group_by(City, ID, country2, Latitude, Longitude) %>%
  filter(Year %in% c(2000, 2019)) %>%
  mutate(change_pm_rate = Rates_PM-lag(Rates_PM),
         change_no2_rate = Rates_NO2-lag(Rates_NO2),
         change_co2 = CO2-lag(CO2),
         City = case_when(City == "Quezon City [Manila]" ~ "Manila",
                          City == "Osaka [Kyoto]" ~ "Osaka",
                          City == "Delhi [New Delhi]" ~ "Delhi",
                          City == "New York" ~ "NYC",
                          City == "Mexico City" ~ "CDMX",
                          City == "Los Angeles" ~ "LA",
                          TRUE ~ City)) %>%
  filter(!is.na(change_pm_rate)) %>%
  dplyr::select(City, ID, country2, Rates_PM, Rates_NO2, CO2, 
                change_pm_rate, change_no2_rate, change_co2)

plot3a <- ggplot(global_cities2, aes(x = reorder(City, -change_pm_rate), y = change_pm_rate, fill = change_pm_rate > 0)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(values = c("steelblue", "salmon"), 
                    labels = c(paste("Decrease"), paste("Increase"))) +
  labs(title = "Top 20 Most Populated Cities by Pct Change in PM2.5",
       x = "City", y = "Change", fill = "Change") 

plot3b <- ggplot(global_cities2, aes(x = reorder(City, -change_no2_rate), y = change_no2_rate, fill = change_no2_rate > 0)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_fill_manual(values = c("steelblue", "salmon"), 
                    labels = c(paste("Decrease"), paste("Increase"))) +
  labs(title = "Top 20 Most Populated Cities by Pct Change in NO2",
       x = "City", y = "Change", fill = "Change") 

pdf(file = "figures/fig3_globalcities.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot3a, plot3b, nrow = 2) 
dev.off()

#-------------------------------------------------------------------------------------------------------------------------------#
# Figure 8: Coal power plant units, installed capacity, and associated emissions added and retired per year for all countries.  #
#-------------------------------------------------------------------------------------------------------------------------------#

dta_plot_coal <- coal_plants %>%
  st_drop_geometry() %>%
  group_by(Start.year) %>%
  summarise(tot_added = n(),
            tot_capacity_added = sum(Capacity..MW.),
            tot_emissions_added = sum(Annual.CO2..million.tonnes...annum.)) %>%
  filter(!is.na(Start.year)) %>%
  filter(Start.year < 2023) %>%
  rename(year = Start.year)

dta_plot_coal2 <- coal_plants %>%
  st_drop_geometry() %>%
  group_by(Retired.year) %>%
  summarise(tot_retired = n(),
            tot_capacity_retired = sum(Capacity..MW.),
            tot_emissions_retired = sum(Annual.CO2..million.tonnes...annum.)) %>%
  filter(!is.na(Retired.year)) %>%
  filter(Retired.year < 2023) %>%
  rename(year = Retired.year)

year <- seq(from = 1927, to = 2022, by = 1)

dta_plot_coal3 <- as.data.frame(year) %>%
  left_join(dta_plot_coal, by = "year") %>%
  left_join(dta_plot_coal2, by = "year") %>%
  mutate_at(vars(tot_added, tot_capacity_added, tot_emissions_added, tot_retired, tot_capacity_retired, tot_emissions_retired), ~replace_na(., 0))

plot5a <- ggplot() +
  geom_line(data = dta_plot_coal3, size = 1.1,
            aes(x = year, y = tot_capacity_added, color = "Added Capacity")) +
  geom_line(data = dta_plot_coal3, size = 1.1,
            aes(x = year, y = tot_capacity_retired, color = "Retired Capacity")) +
  labs(title = "",
       x = "Year",
       y = "Capacity in MW",
       color = "") +
  scale_color_brewer(palette = "Set1") +
  theme_classic()

plot5b <- ggplot() +
  geom_line(data = dta_plot_coal3, size = 1.1,
            aes(x = year, y = tot_added, color = "Added Units")) +
  geom_line(data = dta_plot_coal3, size = 1.1,
            aes(x = year, y = tot_retired, color = "Retired Units")) +
  labs(title = "",
       x = "Year",
       y = "Number of Units",
       color = "") +
  scale_color_brewer(palette = "Set1") +
  theme_classic()

plot5c <- ggplot() +
  geom_line(data = dta_plot_coal3, size = 1.1,
            aes(x = year, y = tot_emissions_added, color = "Added Emissions")) +
  geom_line(data = dta_plot_coal3, size = 1.1,
            aes(x = year, y = tot_emissions_retired, color = "Retired Emissions")) +
  labs(title = "",
       x = "Year",
       y = "CO2 Emissions",
       color = "") +
  scale_color_brewer(palette = "Set1") +
  theme_classic()

pdf(file = "figures/fig5_coalpower.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot5a, plot5b, plot5c, nrow = 3)  
dev.off()

#----------------------------------------------------------------------------------------------------------------------------------#
# Figure 9: Retired coal power plant units, installed capacity, and associated emissions retired per year for selected countries.  #
#----------------------------------------------------------------------------------------------------------------------------------#

dta_plot_coal4 <- coal_plants %>%
  st_drop_geometry() %>%
  group_by(Start.year, Country) %>%
  summarise(tot_added = n(),
            tot_capacity_added = sum(Capacity..MW.),
            tot_emissions_added = sum(Annual.CO2..million.tonnes...annum.)) %>%
  filter(!is.na(Start.year)) %>%
  filter(Start.year < 2023) %>%
  rename(year = Start.year) 

dta_plot_coal5 <- coal_plants %>%
  st_drop_geometry() %>%
  group_by(Retired.year, Country) %>%
  summarise(tot_retired = n(),
            tot_capacity_retired = sum(Capacity..MW.),
            tot_emissions_retired = sum(Annual.CO2..million.tonnes...annum.)) %>%
  filter(!is.na(Retired.year)) %>%
  filter(Retired.year < 2023) %>%
  rename(year = Retired.year) 

coal_country <- data.frame(Country = rep(unique(dta_plot_coal4$Country), each = 99),
                           year = rep(seq(from = 1924, to = 2022), times = 83)) %>%
  left_join(dta_plot_coal4, by = c("Country","year")) %>%
  left_join(dta_plot_coal5, by = c("Country", "year")) %>%
  mutate_at(vars(tot_added, tot_capacity_added, tot_emissions_added, tot_retired, tot_capacity_retired, tot_emissions_retired), ~replace_na(., 0)) %>%
  filter(Country %in% c("United States","China","India","United Kingdom","Indonesia", "Japan", "Germany","Türkiye","Vietnam","South Africa","Vietnam"))

plot6d <- ggplot(coal_country, aes(x = year, y = tot_retired, group = Country, color = Country)) +
  geom_line() +
  xlim(2000,2022) +
  labs(title = "",
       x = "Year",
       y = "Units",
       color = "Country") +
  #scale_color_brewer(palette = "Set1") +
  theme_classic()

plot6e <- ggplot(coal_country, aes(x = year, y = tot_capacity_retired, group = Country, color = Country)) +
  geom_line() +
  xlim(2000,2022) +
  labs(title = "",
       x = "Year",
       y = "Capacity in MW",
       color = "Country") +
  theme_classic() 

plot6f <- ggplot(coal_country, aes(x = year, y = tot_emissions_retired, group = Country, color = Country)) +
  geom_line() +
  xlim(2000,2022) +
  labs(title = "",
       x = "Year",
       y = "CO2 Emissions",
       color = "Country") +
  #scale_color_brewer(palette = "Set1") +
  theme_classic()

pdf(file = "figures/fig6_coalpower_country_retired.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot6d, plot6e, plot6f, nrow = 3)  
dev.off()

#-----------------------------------------------------------------------------------------#
# Figure 10: Total coal-power capacity retired during the period of analysis by country.  #
#-----------------------------------------------------------------------------------------#

dta_plot_coal6 <- dta_plot_coal5 %>%
  group_by(Country) %>%
  summarise(tot_units_retired = sum(tot_retired),
            tot_capacity_retired = sum(tot_capacity_retired)) %>%
  mutate(country2 = countrycode(Country, 'country.name', 'iso3c')) %>%
  filter(!is.na(country2))

plot8b <- ggplot(dta_plot_coal6, aes(x = reorder(country2, -tot_capacity_retired), y = tot_capacity_retired)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  theme_classic() +
  labs(title = "",
       x = "Country", y = "Total Coal Power Retired Capacity (MW)") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1)  # Rotate x-axis text 90 degrees
  )

pdf(file = "figures/fig8_retired_capacity.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot8b
dev.off()

#------------------------------------------------------------------------------------#
# Figure 11: Temporal trends in coal power additions in selected countries.         #
#------------------------------------------------------------------------------------#

plot6a <- ggplot(coal_country, aes(x = year, y = tot_added, group = Country, color = Country)) +
  geom_line() +
  labs(title = "",
       x = "Year",
       y = "Units",
       color = "Country") +
  #scale_color_brewer(palette = "Set1") +
  theme_classic()

plot6b <- ggplot(coal_country, aes(x = year, y = tot_capacity_added, group = Country, color = Country)) +
  geom_line() +
  labs(title = "",
       x = "Year",
       y = "Capacity in MW",
       color = "Country") +
  #scale_color_brewer(palette = "Set1") +
  theme_classic()

plot6c <- ggplot(coal_country, aes(x = year, y = tot_emissions_added, group = Country, color = Country)) +
  geom_line() +
  labs(title = "",
       x = "Year",
       y = "CO2 Emissions",
       color = "Country") +
  #scale_color_brewer(palette = "Set1") +
  theme_classic()

pdf(file = "figures/fig6_coalpower_country.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot6a, plot6b, plot6c, nrow = 3) 
dev.off()

#------------------------------------------------------------------------------------#
# Figure 12: Population living in proximity to a coal power plant.                   #
#------------------------------------------------------------------------------------#

cities_geo_distance2 <- cities_geo_distance %>%
  dplyr::select(1, 4:29)

cities_pop_coal <- dta %>%
  group_by(ID, City, country2) %>%
  summarise(pop = mean(Population, na.rm = T)) %>%
  left_join(cities_geo_distance2, by = "ID") %>%
  mutate(close_operating50k = ifelse(is.na(close_operating50k), 0, close_operating50k)) %>%
  group_by(country2) %>%
  mutate(tot_pop_country = sum(pop)) %>%
  dplyr::select(ID, City, country2, pop, tot_pop_country, close_operating50k) %>%
  ungroup() %>%
  group_by(country2, close_operating50k, tot_pop_country) %>%
  summarise(tot_pop_exposed = sum(pop)) %>%
  mutate(pct_exposed = tot_pop_exposed / tot_pop_country*100) %>%
  filter(close_operating50k ==1) %>%
  filter(pct_exposed > 0)

plot7a <- ggplot(cities_pop_coal, aes(x = reorder(country2, -pct_exposed), y = pct_exposed)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  theme_classic() +
  labs(title = "",
       x = "Country", y = "Percent Population in Proximity") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1)  # Rotate x-axis text 90 degrees
  )

plot7b <- ggplot(cities_pop_coal, aes(x = reorder(country2, -log(tot_pop_exposed)), y = log(tot_pop_exposed))) +
  geom_bar(stat = "identity", fill = "steelblue") +
  theme_classic() +
  labs(title = "",
       x = "Country", y = "Total Population in Proximity (log)") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1)  # Rotate x-axis text 90 degrees
  )

pdf(file = "figures/fig7_pop_exposed.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
grid.arrange(plot7a, plot7b, nrow = 2) 
dev.off()

##############################################
###### Figures Supplementary Information #####
##############################################

#------------------------------------------------------------------------------------#
# Figure S25: Treated and Control Urban Areas United States.                         #
#------------------------------------------------------------------------------------#

xmin <- -66
xmax <- -125
ymin <- 24
ymax <- 49

pdf(file = "si_figures/map1_control_treated.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("USA","Mexico","Canada")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey55",
    alpha = 0.5,
    size = 2.3,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#B10DC9",
    size = 2.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired50k == 1),
    colour = "#FF851B",
    alpha = 0.5,
    size = 2.6,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("operating","construction")),
    colour = "#04ECF0",
    shape = 17,
    size = 1.3,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("retired","mothballed","shelved")),
    colour = "#059DC0",
    shape = 17,
    size = 1.3,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure S26: Treated and Control Urban Areas East Asia.                             #
#------------------------------------------------------------------------------------#

xmin <- 73
xmax <- 145
ymin <- 10
ymax <- 54

pdf(file = "si_figures/map8_control_treated_china.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("China","India","Mongolia",
                                                      "Russia","Thailand","Myanmar",
                                                      "Laos","Japan","Taiwan","Kazakhstan",
                                                      "Vietnam","Kyrgyzstan","South Korea",
                                                      "North Korea","Uzbekistan","Nepal",
                                                      "Bhutan","Pakistan","Cambodia","Philippines")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey55",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#B10DC9",
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired50k == 1),
    colour = "#FF851B",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("operating","construction")),
    colour = "#04ECF0",
    shape = 17,
    size = 1,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("retired","mothballed","shelved")),
    colour = "#059DC0",
    shape = 17,
    size = 1,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure S27: North American Cities by Distance to Retired Coal Power Plants.        #
#------------------------------------------------------------------------------------#

xmin <- -66
xmax <- -125
ymin <- 24
ymax <- 49

pdf(file = "si_figures/map2_distance_retired.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("USA","Mexico","Canada")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey75",
    alpha = 0.5,
    size = 2.3,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired450k == 1),
    colour = "#277da1",
    size = 2.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired150k == 1),
    colour = "#90be6d",
    alpha = 0.5,
    size = 2.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired100k == 1),
    colour = "#f8961e",
    alpha = 0.5,
    size = 2.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired50k == 1),
    colour = "#f94144",
    alpha = 0.5,
    size = 2.6,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("retired","mothballed","shelved")),
    colour = "black",
    shape = 17,
    size = 1.3,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------------------#
# Figure S28: East, South, and Southeast Asian Cities by Distance to Retired Coal Power Plants.  #                          
#------------------------------------------------------------------------------------------------#

xmin <- 73
xmax <- 145
ymin <- 10
ymax <- 54

pdf(file = "si_figures/map9_distance_retired_china.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("China","India","Mongolia",
                                                      "Russia","Thailand","Myanmar",
                                                      "Laos","Japan","Taiwan","Kazakhstan",
                                                      "Vietnam","Kyrgyzstan","South Korea",
                                                      "North Korea","Uzbekistan","Nepal",
                                                      "Bhutan","Pakistan","Cambodia","Philippines")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey55",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired450k == 1),
    colour = "#277da1",
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired150k == 1),
    colour = "#90be6d",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired100k == 1),
    colour = "#f8961e",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired50k == 1),
    colour = "#f94144",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("retired","mothballed","shelved")),
    colour = "black",
    shape = 17,
    size = 1,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure S29: North American Cities by Distance to Operating Coal Power Plants.      #                      
#------------------------------------------------------------------------------------#

xmin <- -66
xmax <- -125
ymin <- 24
ymax <- 49

pdf(file = "si_figures/map3_distance_operating.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("USA","Mexico","Canada")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey75",
    alpha = 0.5,
    size = 2.3,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating450k == 1),
    colour = "#277da1",
    size = 2.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating150k == 1),
    colour = "#90be6d",
    alpha = 0.5,
    size = 2.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating100k == 1),
    colour = "#f8961e",
    alpha = 0.5,
    size = 2.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#f94144",
    alpha = 0.5,
    size = 2.6,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("operating","construction")),
    colour = "black",
    shape = 17,
    size = 1.3,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------------------#
# Figure S30: East, South, and Southeast Asian Cities by Distance to Operating Coal Power Plants   #                         
#------------------------------------------------------------------------------------------------#

xmin <- 73
xmax <- 145
ymin <- 10
ymax <- 54

pdf(file = "si_figures/map11_distance_operating_china.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = subset(world_map, region %in% c("China","India","Mongolia",
                                                      "Russia","Thailand","Myanmar",
                                                      "Laos","Japan","Taiwan","Kazakhstan",
                                                      "Vietnam","Kyrgyzstan","South Korea",
                                                      "North Korea","Uzbekistan","Nepal",
                                                      "Bhutan","Pakistan","Cambodia","Philippines")), aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey55",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating450k == 1),
    colour = "#277da1",
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating150k == 1),
    colour = "#90be6d",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating100k == 1),
    colour = "#f8961e",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#f94144",
    alpha = 0.5,
    size = 1.6,
  ) +
  geom_sf(
    data = subset(coal_plants, Status %in% c("operating","construction")),
    colour = "black",
    shape = 17,
    size = 1,
  ) +
  coord_sf(xlim = c(xmin, xmax), ylim = c(ymin, ymax), expand = FALSE) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure S31: Treated and Control Urban Areas (Global)                               #
#------------------------------------------------------------------------------------#

pdf(file = "si_figures/map6_treated_control_global.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey75",
    alpha = 0.5,
    size = 0.04,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#B10DC9",
    size = 0.04,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired50k == 1),
    colour = "#FF851B",
    alpha = 0.5,
    size = 0.04,
  ) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure S32: World Cities by Distance to Operating Coal Power Plants                #           
#------------------------------------------------------------------------------------#

pdf(file = "si_figures/map5_distance_operating_global.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "grey85", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey65",
    alpha = 0.5,
    size = 0.04,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating450k == 1),
    colour = "#277da1",
    size = 0.04,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating150k == 1),
    colour = "#90be6d",
    alpha = 0.5,
    size = 0.04,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating100k == 1),
    colour = "#f8961e",
    alpha = 0.5,
    size = 0.04,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_operating50k == 1),
    colour = "#f94144",
    alpha = 0.5,
    size = 0.04,
  ) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure S33: World Cities by Distance to Retired Coal Power Plants.                 #           
#------------------------------------------------------------------------------------#

pdf(file = "si_figures/map4_distance_retired_global.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "grey85", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = cities_geo_distance,
    colour = "grey65",
    alpha = 0.5,
    size = 0.06,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired450k == 1),
    colour = "#277da1",
    size = 0.06,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired150k == 1),
    colour = "#90be6d",
    alpha = 0.5,
    size = 0.06,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired100k == 1),
    colour = "#f8961e",
    alpha = 0.5,
    size = 0.06,
  ) +
  geom_sf(
    data = subset(cities_geo_distance, close_retired50k == 1),
    colour = "#f94144",
    alpha = 0.5,
    size = 0.06,
  ) +
  theme(legend.position="none")
dev.off()

#------------------------------------------------------------------------------------#
# Figure S34: Global Fleet of Coal Power Plants.                                     #
#------------------------------------------------------------------------------------#

pdf(file = "si_figures/map7_coalplants.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot() +
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "grey95", color = "black", linewidth = 0.1) +
  theme_minimal() +
  geom_sf(
    data = subset(coal_plants, Status %in% c("operating","construction","retired","mothballed","shelved")),
    colour = "#f94144",
    shape = 17,
    size = 0.6,
  ) 
dev.off()